set.seed(2423045)
rm(list=ls())
library(Matching)

load("rep.sourcecode.re74psid1.RData")

X <-   as.matrix(cbind( foo$age,
                        foo$educ,
                        foo$black,
                        foo$hispan,
                        foo$married,
                        foo$nodegree,
                        I(foo$re74/1000),
                        #I(as.real(foo$re74 == 0)),
                        I(foo$re75/1000)))# ,
                        #I(as.real(foo$re75 == 0))))
                 
BalanceMat <- as.matrix(cbind(
                     foo$age,
                     foo$educ,
                     foo$black,
                     foo$hispan,
                     foo$married,
                     foo$nodegree,
                     I(foo$re74/1000),
                     I(foo$re75/1000),
                     I((foo$re74/1000)^2),
                     I((foo$re75/1000)^2),
                         
                     I((foo$age)^2),
                     I((foo$educ)^2),
                     I((foo$black)^2),
                     I((foo$hispan)^2),
                     I((foo$married)^2),
                     I((foo$nodegree)^2),
                         
                     I(foo$age*foo$educ),
                     I(foo$age*foo$black),
                     I(foo$age*foo$hispan),
                     I(foo$age*foo$married),    
                     I(foo$age*foo$nodegree),
                     I(foo$age*(foo$re74/1000)),
                     I(foo$age*(foo$re75/1000)),

                     I(foo$educ*foo$black),
                     I(foo$educ*foo$hispan),
                     I(foo$educ*foo$married),    
                     I(foo$educ*foo$nodegree),
                     I(foo$educ*(foo$re74/1000)),
                     I(foo$educ*(foo$re75/1000)),

#                     I(foo$black*foo$hispan), This is always zero!
                     I(foo$black*foo$married),
                     I(foo$black*foo$nodegree),
                     I(foo$black*(foo$re74/1000)),
                     I(foo$black*(foo$re75/1000)),

                     I(foo$hispan*foo$married),
                     I(foo$hispan*foo$nodegree),
                     I(foo$hispan*(foo$re74/1000)),
                     I(foo$hispan*(foo$re75/1000)),

                     I(foo$married*foo$nodegree),
                     I(foo$married*(foo$re74/1000)),
                     I(foo$married*(foo$re75/1000)),
                                                
                     I(foo$nodegree*(foo$re74/1000)),
                     I(foo$nodegree*(foo$re75/1000)),                        

                     I((foo$re74/1000)*(foo$re75/1000))                                                 
                         ) )


# From Review of Econ and Statistics
model.A = foo$treat~I(foo$age) + I(foo$age^2) + I(foo$age^3) + I(foo$educ) + I(foo$educ^2) + I(foo$married) + I(foo$nodegree) + I(foo$black) + I(foo$hispan) + I(foo$re74) + I(foo$re75) + I(foo$re74 == 0) + I(foo$re75 == 0) + I(I(foo$re74)*I(foo$educ))

pscores.A <- glm(model.A, family = binomial)

X2 <- cbind(X,BalanceMat[,9:ncol(BalanceMat)]) #first 8 spots of BalanceMat are already in X
for (i in 1:ncol(X2))
  {
    lm1 = lm(X2[,i]~pscores.A$linear.pred)
    X2[,i] = lm1$residual
   }

# VARIABLES WE MATCH ON BEGIN WITH PROPENSITY SCORE
orthoX2.plus.pscore = cbind(pscores.A$linear.pred, X2)
orthoX2.plus.pscore[,1] <- orthoX2.plus.pscore[,1] -  mean(orthoX2.plus.pscore[,1])


# NORMALIZE ALL X COVARS BY STANDARD DEVIATION
#for (i in 1:(ncol(X)+1))
for (i in 1:ncol(X2))
  {
    orthoX2.plus.pscore[,i] <-
      orthoX2.plus.pscore[,i]/sqrt(var(orthoX2.plus.pscore[,i]))
  }



#from psid1.re74.gm4.triton1.Rout
sv <- c(295.511000014901, 14.2662500149012, 0.0082209559011612,
        896.706500014901, 15.0243600149012, 0.228082912914503,
        495.687727940869, 13.4455100149012, 529.198600014901,
        3.59451701490116, 1.49011611938477e-08, 1.30349069030950,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.16125698111802, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 484.152800014901, 1.49011611938477e-08,
        40.8400600149012, 1.49011611938477e-08, 1.75958844981308,
        1.49011611938477e-08, 1.49011611938477e-08, 1.49011611938477e-08,
        1.49011611938477e-08, 1.49011611938477e-08)

#orthoX2.plus.pscore <- orthoX2.plus.pscore[,1:(8+cnew)]

GM.out <-  GenMatch(Tr = foo$treat,
                    X = orthoX2.plus.pscore,
                    starting.values=sv,
                    BalanceMatrix = BalanceMat,
                    pop.size = 1,
                    max.generations = 1,
                    hard.generation.limit=TRUE,
                    wait.generations = 1)

sv=diag(GM.out$Weight.matrix)
cat("Starting Values:\n")
write.csv(sv,eol=", ", row.names=FALSE)

mout <- Match(Y=foo$re78, 
              Tr = foo$treat,
              X = orthoX2.plus.pscore,
              Weight.matrix=GM.out)
summary(mout)


